\documentclass[aps]{revtex4}
\usepackage{graphicx}
\usepackage{amssymb,amsfonts,amsmath,amsthm}
\usepackage{chemarr}
\usepackage{bm}
\usepackage{bbm}
\usepackage{pslatex}
\usepackage{mathptmx}
\usepackage{xfrac}

\newcommand{\mymat}[1]{\bm{#1}}
\newcommand{\mytrn}[1]{~^{\mathsf{T}}{#1}}
\newcommand{\mygrad}{\vec{\nabla}}

\begin{document}
\title{Chemical Solutions}

\section{Description}
Let us assume that we have $A_1,\ldots,A_M$ chemical species coupled by
$N$ equilibria such that
\begin{equation}
	\forall i \in [1;N], \;\; \sum_{j=1}^{M} \nu_{i,j} C_i = 0, \;\; K_i(t) = \prod_{i=1}^{M} C_i^{\nu_{i,j}}.
\end{equation}
We remove the singularities by assuming that the equilibria are met when
\begin{equation}
	\forall i \in [1;N], \;\; \Gamma_i(t,\vec{C}) = K_i(t) \prod_{\nu_{i,j}<0}  C_i^{-\nu_{i,j}} -  \prod_{\nu_{i,j}>0} C_i^{\nu_{i,j}} 
\end{equation}
or
\begin{equation}
	\vec{\Gamma}(t,\vec{C}) = \vec{0}.
\end{equation}
We also naturally have the topology matrix $\mymat{\nu}$, such that
any chemical shift is defined by a  chemical extent vector $\vec{\xi}$ and
produces a shift
$$
	\mytrn{\mymat{\nu}}\vec{\xi}
$$

\section{Limiting extents}
For each equilibrium, depending on the concentration of each components, a 
maximum and a minimum extent may appear, or may
appear unlimited (i.e. for water).
If both way are blocked, then the equilibrium is blocked.



\section{Chemical balancing}
A transformation may bring some active (aka involved in a reaction) concentrations to a negative value.
Let $\vec{C}$ be a set of concentrations. We define the objective function
$$
	\mathcal{E}(\vec{C}) = \sum_j 
	\left
	\lbrace
		\begin{array}{rl}
		-C_j & \text{if species \#}j\text{ is active and } C_j<0\\
		0    & \text{otherwise.}\\
		\end{array} 
	\right.
$$
The steepest integer descent direction is
$$
	\vec{\beta} = -\vec{\nabla}_{\vec{C}} \mathcal{E} = 
	\left\lbrack
		\begin{array}{c}
		\vdots\\
		\mathbbm{1}_{\text{\#}j=\text{active and } C_j<0}\\
		\vdots
		\end{array}
	\right\rbrack
$$
with translate into an integer descent extent
$$
	\vec{\xi}' = \mymat{\nu} \vec{\beta}.
$$
The extent of the blocked equilibria must be set to zero !\\
Then the integer change of concentration direction is
$$
	\delta\vec{C}' = \mytrn{\mymat{\nu}}\vec{\xi}'.
$$
We must then analyse the maximum step that brings us to the
first zero concentration, while taking care of the addition and subtraction rounding.
If this step is successful, then at least one of the negative species is exactly zero...\\

The use of this $\mathcal{E}$ make a non vanishing $\delta\vec{C}'$.

\section{Newton's algorithm}
Let us start from a balanced concentration $\vec{C}_k$.
We want to find
$$
	\vec{\Gamma}(\vec{C}_k + \mytrn{\mymat{\nu}}\vec{\xi}_k) = \vec{0}.
$$
The first order expansion provides
$$
	\vec{\Gamma}(\vec{C}_k + \mytrn{\mymat{\nu}}\vec{\xi}_k) \simeq 
	\vec{\Gamma}_k + \mymat{\Phi}_k\mytrn{\mymat{\nu}} \vec{\xi}_k
$$
For a valid composition, the Newton's extent is
$$
		\vec{\xi}_k = -\left(\mymat{\Phi}_k\mytrn{\mymat{\nu}}\right)^{-1}\,\vec{\Gamma}_K
$$
which must be clamped to deduce the modified Newton's step
$$
	\delta{\vec{C_k}} = \mytrn{\mymat{\nu}} \vec{\xi}_k^\star
$$
and a careful addition is then performed, avoiding the roundoff errors.

\section{Booting Algorithm}
\subsection{The problem}
We have $M$ species, $N$ equilibria, 
so we must have $M\geq N$ and $Nc=M-N$ other conditions
to obtain a starting value.
Let us assume that we deal only with linear constraints (electroneutrality, mass conservation...).
This may be written in terms of an integer matrix $\mymat{P}$
such that
$$
	\mymat{P}\vec{C} = \vec{\Lambda}.
$$
To be able to find a solution, the matrix $\mymat{P}\oplus\mymat{\nu}$ must
be not singular, with allows the computation of an orthogonal (integer) matrix perpendicular to $\mymat{P}$
$$
	\mymat{Q},\;\mytrn{\mymat{P}}\cdot\mymat{Q} = \mymat{0},\;Q\in\mathcal{M}_{N,M}.
$$
Any valid solution must be of the form
$$
	\vec{C}_k = \vec{C}^\star + \mytrn{\mymat{Q}}\vec{\xi}_k
$$
with
$$
	\vec{C}^\star = \mytrn{\mymat{P}}\left(\mymat{P}\mytrn{\mymat{P}}\right)^{-1} \vec{\Lambda}
$$
\subsection{Best effort positive values}

\section{Perturbation}
Let us assume we want to add $\delta\vec{C}$ to the system during $dt$, starting
from $\vec{C}$ with $\vec{\Gamma}(\vec{C})=\vec{0}$.
We want to compute the extent $\delta\vec{\xi}$ such that
$$
	\vec{\Gamma}\left(\vec{C}+\delta\vec{C} + \mytrn{\mymat{\nu}}\delta\vec{\xi}\right)=\vec{0}
$$
imposing the first order relation
$$
	\mymat{\Phi}\delta\vec{C} + \mymat{\Phi}\mytrn{\mymat{\nu}} \delta\vec{\xi}
	+\partial_t \vec{\Gamma} dt = \vec{0}.
$$
using 
$$
	\vec{\rho} = \partial_t \delta\vec{C}
$$
we find
$$
	\dot{\vec{\xi}} = - \left(\mymat{\Phi}\mytrn{\mymat{\nu}}\right)^{-1}
	\left(
		\mymat{\Phi}\vec{\rho} + \partial_t \vec{\Gamma}
	\right)
$$
so that 
$$
	\delta\vec{\rho} = - \mytrn{\mymat{\nu}} \left(\mymat{\Phi}\mytrn{\mymat{\nu}}\right)^{-1}
	\left(
		\mymat{\Phi}\vec{\rho} + \partial_t \vec{\Gamma}
	\right)
$$

\end{document}